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SUMMARY 

This  is  a written  version  of  a lecture  given  by  the  author  at  Euromech 
Colloquium  75  which  was  held  at  Rhode,  near  Braunschweig  in  May  1976.  The 
subject  of  the  colloquium  was  'The  calculation  of  flow  fields  by  panel  methods' 
and  this  paper  describes  iterative  calculations  of  the  classic  problems  of 
steady  inviscid  flow  around  a thick  cambered  wing,  providing  improved  accuracy  in 
relation  to  the  so-called  RAE  Standard  Method. 
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1 INTRODUCTION 

Over  a period  of  many  years,  an  approximate  method  for  the  solution  of  the 
classical  problem  of  steady  inviscid  flow  around  a thick  cambered  isolated  wing 
was  developed  in  England,  principally  by  Dr  Weber  and  the  late  Dr  Kuchemann,  and 
became  known  as  the  RAE  Standard  Method.  This  method  is  based  effectively  on 
the  representation  of  the  wing  by  planar  source  and  doublet  distributions;  in 
first-order  theory,  sources  represent  wing  thickness  and  doublets  represent  wing 
loading.  The  strengths  of  these  planar  singularity  distributions  could  be 
estimated  by  appeal  to  the  theory  of  infinite  swept  wings  and  of  infinite  kinked 
wings,  with  some  interpolation  between  relevant  regions.  This  method  was  rela- 
tively fast,  and  as  it  was  based  on  the  physics  of  the  problem,  it  was  relatively 
easy  to  understand.  However,  in  time  questions  were  raised  whether  the  estimates 
involved  were  sufficiently  accurate.  One  way  to  answer  these  questions  was,  to 
write  some  computer  programs  to  evaluate  accurately  the  velocity  fields  due  to 
given  planar  singularity  distributions,  by  performing  the  necessary  double  inte- 
grations, and  in  1968  to  1970  such  programs  were  developed  at  RAE  by  John  Ledger ^ 

2 . . 
and  the  author  . Separate  versions  were  written  to  calculate  the  flow  fields  on 

and  off  the  singularity  plane.  Next,  in  1972  Johanna  Weber  showed  how  to  apply 

the  results  of  these  computations  iteratively  to  improve  the  basic  solution  to  at 

least  second-order  accuracy.  In  the  following  years  the  author  wrote  a further 
4 5 

set  of  programs  ’ to  implement  her  ideas,  and  it  is  these  iterative  calculation 
methods  that  are  described  here.  The  whole  group  of  methods  thus  represents  a 
refinement  of  the  RAE  Standard  Method,  although  because  of  the  numerical  double 
integrations  there  is  a price  to  be  paid  in  computing  time. 

2 DIRECT  PROBLEMS 

Fig  1 shows  half  of  a symmetrical  wing  at  incidence  a . We  define  the 
cartesian  coordinate  y measured  to  starboard  from  the  plane  of  symmetry,  and 
at  each  y we  consider  a plane  containing  the  local  section  chordline  and  at 
incidence  a ♦ ctj,  to  the  free  stream,  (*T  being  the  local  twist.  The  singular- 
ity distributions,  which  are  actually  in  the  wing  chordal  surface,  can  be  thought 
of  as  lying  in  this  plane,  to  second-order  accuracy.  We  can  take  local  cartesian 
coordinates  (x,z)  with  the  x-axis  parallel  to  the  local  chordline  and  the  z-axis 
upwards,  completing  the  right-handed  set.  Wing  surface  ordinates  in  this  system 
now  vanish  at  the  leading  and  trailing  edges.  We  define  them  thus: 


zw(x»y)  ■ ± zt(x,y)  + zg(x,y)  • 


z * 
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z^  is  the  thickness,  zg  the  camber  ordinate;  upper  and  lower  signs  correspond 
to  upper  and  lower  wing  surfaces.  For  the  local  incidence  a + , the  free 

stream  flow  is  normalized  to  unit  speed  and  written: 


U 


[cos  (a  + c^),  0,  sin  (a  + c^)] 


In  the  direct  problem4,  we  aim  to  satisfy  the  boundary  condition  of  zero 
normal  flow  by  iterative  calculation  of  source  and  doublet  distributions  q and 
l . Let  us  suppose  that  we  have  obtained  a set  of  trial  singularity  distribu- 
tions. We  now  calculate  their  separate  velocity  fields  using  the  computer  sub- 
routines written  by  Ledger  and  the  author. 

For  an  uncambered  wing,  z = 0 , we  can  calculate  the  velocities  on  the 

s 

upper  surface  z = z^  and  use  these  relations  for  the  source  field  ut  and 

doublet  field  u„ 

-l 


*t 


-H 


(ut»  vt,  +wt) 


(±V  ±V£*  V • 


to  obtain  the  values  on  the  lower  surface.  In  general,  for  a cambered  wing,  we 

still  compute  velocity  fields  on  the  thickness  surface  z = z^  and  obtain  the 

values  on  the  actual  wing  surface  by  Taylor  series  expansions  in  z : for 

s 

example,  the  streamwash  ufc  due  to  sources  is  written: 


3u 

ut(x,y,zt±zs)  * ut(x,y,zt)  ± zs  (x,y,zt) 

with  similar  expressions  for  the  other  five  components. 

Substituting  all  these  into  the  boundary  condition 


R = + “t  + • grad  ( z w - z) 


we  get 
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We  multiply  this  out  and  take  the  plus  or  minus  parts  to  Rt  and  the  rest  to 

R . In  the  special  case  of  an  uncambered  wing,  z * 0 , assuming  small  inci- 
X.  s 

dence,  these  expressions  simplify  to 


3zt  3zt 

(1  +ut)  3ir  + vt  37" -wt 


3zt  3zt 

u*  TT  * aT  " w*  ' (a  + aT) 


In  this  simple  case,  the  source  and  doublet  problems  uncouple.  In  general, 
however,  the  problems  must  be  solved  together. 

These  residual  errors  are  evaluated,  and  tell  us  not  only  how  well  the 
boundary  conditions  have  been  satisfied,  but  also  what  perturbation  source  and 
doublet  distributions  to  add  on  in  order  to  cancel  R^  and  R^  approximately. 
Rt  is  regarded  as  a shortfall  in  wt  (if  Rt  is  positive,  w£  is  not  big 
enough),  and  similarly  R is  regarded  as  a shortfall  in  w . So  to  cancel 

X-  X, 

Rt  we  take 

Aq  - 2Rt  , 

and  we  generate  the  extra  doublet  distribution  from  R^  by  an  approximate  but 
rapid  lifting-surface  theory.  We  use  a vortex-lattice  method  for  that  job.  The 
iteration  cycle  is  now  complete,  and  can  be  repeated  as  desired. 


Starting  values  of  q and  l can  be  obtained  from  linear  theory: 
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For  subcritical  flow  at  low  Mach  number  M , we  can  transform  to  affine 

00 

or  Prandtl-Glauert  variables: 


xB 


ut/B 


u*  “ VB 


(B2  - 1 - M2) 


Then  in  the  affine  (x,y,z)  space  we  again  have  an  incompressible  flow  which  can 
be  built  up  from  source  and  doublet  fields.  Note  that  this  accounts  for  linear 
compressibility  effects  only.  Note  also  that  the  flow  is  not  exactly  equivalent 
to  flow  over  an  analogous  wing,  because  the  boundary  conditions  are  slightly 
different.  For  example,  considering  the  uncambered  wing  again,  is  given 

thus : 

/ \\  l 3zt  3zt 

Rt  y + B J B 3x  + Vt  3y  ~ Wt 

There  are  two  extra  factors  B present. 

The  computation  of  velocity  fields  on  the  wing  surface,  even  at  a finite 
number  of  collocation  points,  is  lengthy.  So  we  would  like  to  keep  the  number 
of  iteration  cycles  as  small  as  possible.  One  trick  is  to  expand  the  residual 
errors  R and  R as  Maclaurin  series  (about  z **  0)  in  z thus: 

t X.  W 


(n,p+l) 


cos  (a  + a„)  + u(x,y,z  ) + Au(x, 
1 w 


]9Z 

sr 

1 9z 

’*‘Zw>J  3F 

- |sin  (a  ♦ aT)  + w(x,y,*M)  + Aw(x,y,zw)J 
8<».P>  . (to  . z.  |f)  (av  . z„  ff)  £ - (A„  . z„  |f) 


v(x,y,z  ) + Av(x, 
w 


All  quantities  are  now  evaluated  on  z » 0 . p is  an  inner  iteration  super- 
script. Since  (in  incompressible  flow) 


9 Aw  ^ _ 3Au  3Av 

9z  3x  9y 


,(n,p+l ) 


.(n.P) 


- Aw  + 


97  (ZWAU) 


t—  (z  Av) 
9y  w 


we  have 
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Substituting  everything  in,  multiplying  out  and  separating  the  ± parts,  we  get 
these  two  expressions: 


7 


rj 


R 


(n,p+l)  _ R(n,p)  - 


lut  * £ (ztaut  * W * F (V\  ♦ W 


R 


(n,p+l)  _ n(niP)  _ 


Aw. 


+ £ 


z Au  ) + r — (z  Av„ 
st  iy  ' t t 


+ z Av_) 
s t 


The  first  two  terms  on  the  right  sides  have  been  cancelled  by  the  iterative 

selection  of  source  and  doublet  distributions.  We  can  now  use  approximate 

expressions  for  Au  , Av  in  terms  of  Aq  and  for  Au  , Av.  in  terms  of  AJ,  , 
t t a a 

on  z = 0 , instead  of  computing  them  accurately. 

For  simple  wings,  two  main  iterations  often  suffice.  For  curved  wings,  or 
wings  with  cranks,  or  complex  camber  shapes,  three  or  more  iterations  should  be 
run. 

Some  comparisons  have  been  made  with  the  BAC  Neumann  program  of  Roberts, 

the  standard  datum  method,  for  the  first  family  of  uncambered  wings  taken  by 

7 8 

BAC  and  NLR  as  a standard  test  case,  and  discussed  by  Hunt  (BAC)  ’ . These  wings 
have  the  planform  of  RAE  Wing  'A',  aspect  ratio  6,  taper  ratio  3,  mid-chord  sweep 
angle  30°,  and  the  chordwise  thickness  distribution  is  that  of  the  NACA  00  series. 
Fig  2 shows  some  results  for  a wing  with  a 15%  thick  section,  at  zero  incidence, 
at  three  stations  across  the  semispan.  Agreement  is  good,  except  perhaps  near 
the  thick  trailing-edge ; Hunt  has,  however,  now  told  us  that  Roberts  had 
extended  this  to  a point.  The  top  part  of  Fig  3 shows  some  results,  for  the  same 
wing  at  5°  incidence;  there  is  a slight  discrepancy  near  the  leading-edge  in  mid- 
semispan. Perhaps  the  author's  method  needs  more  than  11  collocation  points.  On 
the  bottom  part  of  Fig  3 are  results  for  the  much  thinner  wing,  2%  thick,  at  the 
same  incidence  and  spanwise  station.  There  is  no  leading-edge  discrepancy  here. 


Comparisons  with  the  Roberts  method  have  also  been  made  for  RAE  Wing  'B', 
of  which  the  planform  (the  same  as  that  of  Wing  'A'),  camber  and  twist  distribu- 
tions are  shown  in  Fig  A.  For  zero  incidence  and  zero  Mach  number,  the  results 

in  mid-semispan  compare  well  (Fig  5) . We  also  show  results  from  a fast  approxi- 

. 4 

mate  method  due  to  Lock.  Comparisons  are  not  so  good  near  the  root,  where 
Wing  'B'  has  local  dihedral  which  the  present  method  neglects,  and  near  the  tip 
where  neither  the  present  method  nor  Roberts'  method  satisfies  the  end  boundary 
condition  on  the  square  cut  tip. 


k.  k 
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Ground  effect  can  be  represented  by  image  source  and  doublet  distributions, 
of  the  same  sign  for  images,  of  opposite  sign  for  doublets.  The  iteration 
proceeds  much  as  before;  the  velocity  fields  due  to  the  image  distributions  are 
also  evaluated  by  the  Ledger-Sells  subroutine  on  the  real  wing  chordal  surface, 
and  combined  as 


- (uG‘VG*WG)  • 

The  boundary  condition  now  reads 


R'  = <!L  + yt  + “a  + yG)  • grad  (z^  - z)  - 0 . 

As  before,  this  can  be  split  as  ± Rfc  + R^  = 0 . For  zfi  = 0 again,  we 

Rt  = (1  + \ + ug)  jt  + (vt + v 3 r ~ wt 


3z  3 z 

K = u*  air + \ ar  - w*  - wG  - (a  + V • 


For  simple  wings  the  program  now  needs  three  iterations  instead  of  two,  probably 
because  the  perturbation  singularity  distributions  are  generated,  neglecting  the 
effect  of  their  own  images. 

In  the  current  version,  we  have  neglected  variations  in  u between  upper 

*“G 

and  lower  surfaces  for  each  (x,y) . Comparison  with  results  from  the  Hess  and 
Smith  surface  panel  method  for  a two-dimensional  wing  suggests  that  this  may  have 
been  unwise  . 

3 DESIGN  PROBLEMS 

Programs  based  on  this  method  have  also  been  written  for  various  design 
options'*.  In  order  of  description,  we  can  specify  the  thickness  and  the  first- 
order  loading  (that  is  to  say,  the  doublet)  distributions;  or  we  can  specify  the 
thickness  and  the  upper-surface  pressure  distribution;  or  we  can  specify  the 
first-order  loading  and  the  upper-surface  pressure  distribution.  In  each  case, 
the  camber  and  twist  distributions  are  determined  as  part  of  the  solution. 

In  all  cases,  the  source  distribution  has  to  be  found,  and  this  is  done 
iteratively  using  the  successive  error  field  iterates  R^  exactly  as  before. 
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The  successive  iterates  for  the  other  error  field 


are  now  used  to 


generate  corrections  to  the  camber  and  twist  distributions  A z , Aa  . If  we  go 

back  to  the  boundary  condition,  substitute  the  new  camber  and  twist  distributions 

z + Az  , a + Aa„  , neglect  products  of  perturbation  quantities  and  equate  the 
S S 1 X 

resulting  expression  to  zero,  we  get: 


Since  wing  ordinates  vanish  at  leading  and  trailing  edges  in  our  coordinate 
system,  integrating  over  the  whole  chord  gives  the  wing  twist: 

*T  / 

Ac*T(y)  = j R^Cx’  ,y)dx7  [x^y)  - x^y)]  . 

*L  / 

Then  the  indefinite  integral  gives  the  camber  correction: 

x 

Azg(x,y)  = [x  - xL(y)]AaT(y)  - J R^(x’,y)dx’  . 


This  is  all  that  is  needed  for  Option  1,  in  which  thickness  and  first- 
order  loading  are  specified.  The  upper-surface  pressures  come  out  as  part  of 
the  solution.  This  option  converges  quickly,  about  as  quickly  as  the  direct 
program. 

For  Option  2,  when  we  specify  thickness  and  upper-surface  pressure  , 

we  now  have  to  calculate  the  doublet  strength  to  satisfy  the  condition 

iteratively.  A first  estimate  for  £,  is  provided  by  a modification  of  the  RAE 
Standard  Method  due  to  Lock.  In  subsequent  iterations,  the  shortfall  in  C 

pu 

is  in  principle  a second-order  quantity  and  from  it  we  can  derive  a second-order 
expression  for  the  correction  to  the  streamwash  Au^  due  to  loading.  We  will 
use  suffix  u to  represent  upper-surface  values.  The  current  local  speed  Q 
is  given  by 

2 2 2 2 

Q = u + V + W . 
u u u u 

For  the  design  speed  we  can  write: 


^ It 
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-2  2 2 2 
Q = (U  + Au.)  + V + W 
u V u u 


Q + 2U  Au„ 
u u l 


From  this  equation  Au^  follows: 


Au  = (Q2  - Q2)/2U 

i u u j u 


Then  Ai.  is  given  thus: 


A£  = 46Au^ 

This  iteration  cycle  converges  in  mid-setnispan,  and  with  the  help  of  some 
relaxation  factors  it  also  converges  near  the  tip,  for  the  cases  studied.  But 
near  the  root  of  a swept  wing  convergence  is  slow,  and  the  iteration  may  in  fact 
diverge.  It  seems  that  near  the  root  the  logarithmic  singularity  of  linear 
theory  in  the  upwash  Aw^  affects  the  velocity  components  so  that  a small  change 
Ail  does  not  necessarily  produce  small  changes  in  v^,  w^  and  so  does  not  have 
the  desired  effect  on  . For  some  cases,  we  avoided  instability  by  arranging 

an  inner  iteration  cycle  in  which  we  calculate  Av^  and  Aw^  approximately  and 
then  update  and  generate  a further  correction  doublet  field  as  before. 

This  also  speeds  up  convergence  over  the  rest  of  the  wing. 


The  program  was  applied  to  RAE  Wing  ' B'  at  M = 0.8  , as  a check.  Target 


distributions  of  H , C 


pu 


section  lift  C and  centre  of  pressure  x were 
LL  cp 


generated  using  the  earlier  design  program  (for  Option  1),  and  are  shown  as  full 
lines  in  Fig  6.  The  top  part  shows  the  convergence  near  the  root  of  the  load 
function  (this  is  a function  made  regular  at  the  leading  and  trailing  edges  by 
multiplying  out  the  square  root  singularities).  The  load  function  has  not  quite 
converged  near  the  root.  However,  in  mid-semispan  and  near  the  tip  convergence 
is  excellent^.  This  is  reflected  in  the  lower  parts  of  Fig  6,  where  the  twist 
and  the  section  aerodynamic  quantities  are  very  well  converged  outboard,  but  not 
quite  on  target  near  the  root.  Fig  7 tells  the  same  story:  the  camber  distribu- 
tion is  well  converged  near  the  tip  and  in  mid-semispan,  but  has  some  way  to  go 

near  the  root.  Fig  8 shows  the  same  again  for  C , but  we  see  that  the  final 

pu 

iterates  are  very  close  to  the  target  near  the  root,  in  marked  contrast  to  the 
camber  and  twist,  so  that  C near  the  root  is  evidently  not  very  sensitive  to 


the  wing  shape, 


pu 


Here, 


For  Option  3,  we  specify  upper-surface  pressure  and  first-order  loading. 

it  is  the  thickness  that  must  be  adjusted  to  get  the  right  C 

t pu 
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I 


L. 


distribution.  This  time  the  shortfall  in 
velocity  due  to  sources  thus: 


C 

pu 


is 


-2  2 

Au  = (Q  - (H/2U 
t u xu  u 


converted  into  a shortfall  in 


The  right  side  of  this  equation  is  the  same  as  that  for  the  doublets  in  Option  2. 
The  required  perturbation  thickness  Az^  can  be  thought  of  as  giving  an  extra 
source  distribution  28Azt/9x  . From  the  RAE  Standard  Method  we  have  the 
following  approximate  formula  for  Au^  in  terms  of  Azt  (the  simplified  rela- 
tion for  incompressible  flow  is  given  here) : 


Au„  = cos  A 


1 

n 3x' 

j 


9AZt(x',y)  dx, 


x - X 


v (y)  I In  .t-4n  K- 
2 J it  1 - sin  A 


3Az{ 

3x 


with  A = local  sweep  angle 

= spanwise  interpolation  parameter. 


1 + 


'9Azt 

9x 


In  a sheared-wing  region,  or  near  mid-semispan,  K2  = 0 and  this  is  an 
ordinary  Cauchy-type  principal-value  integral  equation,  and  the  solution  is 
standard.  When  K2  is  not  zero,  we  can  solve  this  equation  iteratively,  by 
regarding  the  last  term  as  known  from  a previous  iterate,  and  again  inverting 
the  standard  integral  equation. 


This  option  converges  slowly  near  the  root,  but  reasonably  well  elsewhere. 

It  is  possible  to  combine  the  second  and  third  options  into  a hybrid 
Option  4.  In  this  option,  a station  y = y*  is  picked  near  the  root.  Outboard 
of  y = y*,  C and  z t are  specified,  leading  to  a sequence  of  iterates  for 
the  outboard  doublet  distribution  as  in  Option  2,  and  then  the  doublet  distribu- 
tion is  constrained  to  exhibit  the  same  ch^rdwise  behaviour  throughout  the  region 

inboard  of  y = y*  , up  to  the  root.  To  obtain  a specified  C in  this  region 

pu 

also,  we  can  now  adjust  the  inboard  thickness  distribution  z^  as  in  Option  3. 
This  option  thus  avoids  the  convergence  problems  associated  with  doublets  near 
the  root. 

As  a test  case  for  this  option,  a modified  wing  was  first  designed  to 

have  a suitable  load  distribution  and  RAE  101  thickness  distribution,  9%  thick 

outboard,  rising  to  13.5%  at  the  root.  This  gave  an  output  C distribution 

pu 
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dis tributions  of  z_/c  near  the  root,  of  wing  twist  and  of  CTT  and  x . 

t ll  cp 

The  first  guess  at  z^/c  was  the  same  as  that  outboard,  the  91  thick  RAE  101 

distribution.  After  four  more  iterations  the  inboard  thickness  is  slowly 

converging,  but  is  not  quite  on  target  yet,  and  we  expect  this  to  show  up  when 

comparing  other  quantities.  Indeed,  in  the  lower  part  of  the  figure  the  twist 

and  the  section  lift  and  centre  of  pressure  are  well  converged  outboard,  but  not 

quite  on  target  near  the  root.  Fig  10  shows  the  camber  iterates:  again  well 

converged  outboard,  but  not  quite  home  near  the  root.  Fig  11  shows  the  upper- 

surface  pressures:  also  well  converged  outboard,  and  slight  differences  near  the 

root.  In  fact,  C does  not  seem  very  sensitive  to  root  thickness, 
pu 

It  can  be  seen  that  this  Option  4 gives  the  designer  a chance  to  maintain 
good  flow  quality  right  into  the  root  of  a swept  wing,  by  increasing  the  thick- 
ness which  he  wants  to  increase  anyway  to  strengthen  the  wing-body  junction. 
Thus,  it  seems  that  this  could  be  a quite  useful  design  option. 
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Fig. I Sketch  of  half  wing  and  typical  singularity  plane 
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Fig.  5 Comparison  of  methods  at  mid-semi-span  for  RAE 

wing'B'  at  zero  incidence 


Target 

— 

Iteration  1 

+ (±) 

Iteration  4 

© 

T.Motmo  Aero  1709 


r.  Memo  Aero  1703 


0050r 


iL 

c 


0 0-2  0-4  0-6  0-8 

a Thickness  distribution  near  root 


Fig.  9 


Target  (known) 

— 

Starting  distribution  + 

Iteration  2 

X 

1 teration  5 

o 

0 0-2  0-4  0-6  0 8 1-0 

c Section  lift  and  centre  of  pressure 


Fig. 9 Iterates  for  RAE  wing  “B" Moo  = 0*8.  Problem  4 


Fig. 10  Iterates  for  camber  distribution  on 
R AE  wing  "B",  = 0 • 8.  Problem  2 
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Fig  II  Iterates  for  upper-surface  pressure  distribution  on 
R AE  wing  " B",  Mod  = 0*8.  Problem  4 
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